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ABSTRACT 


A method, called the Eigensystem Realization Algorithm ( EF 5 developed 
for modal parameter identification and model reduction of dynamic systems from 
test data. A new approach is introduced in conjunction with the singular value 
decomposition technique to derive the basic formulation of minimum order 
realization which is an extended version of the Ho-Kalman algorithm. The basic 
^ormulation is then transformed into modal space for modal parameter Identifica- 
tion. Two accuracy indicators are developed to quantitatively Identify the 
system modes and noise modes. For illustration of the algorithm, examples are 
shown using simulation data and experimentaldata for a rectangular grid struc- 
ture. 


I. INTRODUCTION 

The state space model has received considerable attention for system 
analyses and design In recent control and systems research programs One of 
these areas. In particular. Is control of large space structures. In order to 
design controls for a dynamic system It Is necessary to have a mathematical 
model which will adequately describe the system's motion. The process of 
constructing a state space representation from experimental data Is called 
system realization. 

During the past two decades, numerous algorithms for the construction of 
state space representations of linear systems have appeared In the control 
literature. Among the first were the works of Gilbert tl] and Kalman [2], 
introducing the important principles of realization theory In terms of the 
concepts of controllability and observability. Both techniques use the transfer 
function matrix to solve the realization problem. Ho and Kalman [3] approached 
this problem from a new viewpoint. They showed that the minimum realization 
problem is equivalent to a representation problem Involving a sequence of real 
matrices known as Markov parameters (pulse response functions). By nlnlmum 
realization is meant the smallest state-space dimension among systems realized 
that has the same input-output relations within a specified degree of accuracy. 
Questions regarding the minimum realization from various types of Input-output 
data and the generation of minimum partial realization are studied by Tether 
[4], Silverman [5], and Rossen and Lapldus [ 6 ] using Markov ^■'’-neters. Rossen 
and Lapldus [7] successfully applied Ho-Kalman [3] and Tether [4] methods to 
chemical engineering systems. A common weakness of the above schemes Is that 
effects of noise on the data analysis were not evaluated. Zelger and McEwen 
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[8] proposed a combinition of the Ho-Kalman algorithm [3] with the singular 
value decomposition technique for the treatment of noisy date However, no 
theoretical or numerical studies were reported in Reference [8 1 vnong follow- 
up developments along similar lines, Kung [9] presented anc algorithm in 
conjunction with the singular value decomposition technique V'. /ncorporate the 
presence of the noise. Note that the singular value decoir.x^ition technique 
[10-11] has been widely recognized as being very effective and numerically 
stable. Although several techniques of minimum realization are available in the 
literature, formal direct application to the modal parameter identification for 
flexible structures was not yet addressed. 

In the structures fieiu, the finite-element technique is used almost 
exclusively for constructing analytical models. This approach is well estab- 
lished and normally provides a model accurate enough for structural design 
purposes. Once the structure is built, static and dynamic tests are performed. 
These tes. results are used to refine the finite-element model, which is then 
use for-afinal analyses. This traditional approach to analytical model devel- 
opment may not be accurate enough for use in designing a vibration control 
system for flexible structures. Another approach is to realize a model directly 
from the experimental results. This requires the construction of a minimum- 
order model from the test data that characterizes the dynamics of the system 
at the selected control and measurement positions. The present state-of-the-art 
in structural modal testing and data analysis is one of controversy about the 
best technique to use. Classical test techniques, which may provide only good 
frequency and moderate mode shape accuracy, are often considered adequate for 
finite-element model verification purposes. On the other hand, advanced data 
analysis techniques which offer significant reductions in test time and Improved 
accuracy , have been available [12-16] but are not yet fully used, for example, 
Ibrahim [13] presented a method based on state space for the direct identifica- 
tion of modal parameters from free responses. Recently, Void and Russell [16] 
presented a method using frequency response functions and time domain analysis 
which can also identify repeated eigenvalues. A comparison of contemporary 
methods using data from tne Galileo spacecraft test Is provided by Chen [17]. 

Although structural dynamics techniques are generally successful for ground 
data, further incorporation with work from t.he controls discipline is needed to 
solve modal parameter identification/control problem. For example, it is known 
from control theory [18] that a system with repeated eigenvalues and Independent 
mode shapes is not identifiable by single input and single ouput. Methods 
which allows only one Initial condition (Input) at a time [13], will miss 
repeated eigenvalues. Also, if the realized system is not of a minimum order 
and matrix inversion is used for constructing an oversized state matrix, 
numerical errors may become dominant. 

Under the interaction of structure and control disciplines, the objective 
of this paper is to introduce an Eigensyscem Realization Algorithm ( ERA ) for 
modal parameter identifica <on and model reduction for dynanneal systems from 
test data. The algorithm consists of two major parts, namely, basic formulation 
of the minimum order realization and modal parameter identification. In the 
section of basic formulation, the Hankel matrix which represents the data 
structure for Ho-Kalman algorithm is generalized to allow random distribution 
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of Markov parameters generated by free decay reponses. A unique approach 
based on this generalized Hankel matrix is developed to extend the Ho-Kalman 
algorithm in combination with the singular value decomposition technique [10- 
11] . Through the use of the generalized Hankel matrix, a linear model is realized 
for dynamical system matching the input and output relationship. The realized 
system model is then transformed into modal space for modal parameter identifi- 
cations. As part of ERA, two accuracy indicators, namely, the modal amplitude 
coherence and the modal phase col linearity, are developed to quantify the 
system modes and noise modes. The degree of modal excitation and observation 
are evaluated. The ERA method thus forms the basis for a rational choice of 
model size determined by the singular values and accuracy indicators. 

Two examples are given to illustrate the ERA method. The first example 
uses simulated data from an assured structure. The effect of repeated eigenvalues 
on the parameter identification is shown. The second example uses experimental 
data from \ simple grid structure. Comparison of the ERA results with a finite 
element model of the grid is performed. Experimental results for a more complex 
structure— the Galileo spacecraft—are shown in Ref. [19]. 


II. BASIC FORMULATIONS 


A finite-dimensional, discrete time, linear, time-invariant dynamical 
system has the state-variable equations 


x(k+l) * Ax(k) + Bu(k) 

(1) 

y(k, = Cx(k) 

(2) 


where x is an n dimensional state vector, u is an m dimensional control input, 
and y is an p dimensional output or measurement vector. The integer k is the 
sample indicator. The transition matrix A characterizes the dynamics of the 
system. For flexible structures, it is a representation of mass, stiffness 
and damping properties. The problem of system realization is then the follow- 
ing. Given the measurement functions y(k), construct constant matrices [A, B, 
C] such that the functions y are reproduced by the state-variable equations. 
With different sets of inputs and outputs, several cases can be obtained. The 
simplest case, namely, single input and single output, is treated first to 
allow the reader familiar with notations for the treatment of multi-input and 
multi -output cases. 

Single input and single output 

For the system (1) with free pulse-response (or initial-state-response), 
the time-domain description is given by the function known as Markov parameter 

y(k) - CA k "lB [or y(k) » CAkx(O)] (3) 

where x(0) is the system ini tie' onditions and k is an integer. Note that the 
matrix B is a column vector (si .gie input) whereas the matrix C is a row vector 
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(single output). For free init1al-state>respon$e, the matrix B only represents 
the information of initial conditions rather than the control influence matrix 
as shown in Eq.(l). The problem of system realization is to construct matrices 
[A, B, C] in terms of the measurement function y(k) such that the identities of 
Eq.(3) hold. Now observe that 


where 


y(k) * VA k- lB [ or y(k) « VA k x(0) ] 
y(k) 



ryOO 1 


C ' 

* 

y(k+i) 

• 

; and V = 

CA 

• 


• 

_y(k+n-l)_ 


• 

CAd-1 


(4) 

(5) 


Assume that this nth order system has no repeated eigenvalues. There exists 
a row vector C from observability theory (Ref. 18) such that V has rank n. 
Consequently, rearranging Eq. (4) becomes 

y(k+l ) = VA*B * VAV^y(k) (6) 


Given the sequence of measurement vectors y(k+l), 
matrix H(k) is defined as 


H(k-l) * Cy(k) t y(k+l) 


y(k+n)] 


fy(k) M™) 

y(k+l) ,y(k+2) 


the general i zed 

..... y(k+n-l)l 
..... y(k+n) 


Hankel 


(7) 


[y(k+n-l) ,y(k+n) y(k+2n-2)j 


It immediately follows from Eq.(6) that 

H(k) = VAV'^K(k-l) » VA k V“ 1 M(0) 


( 8 ) 


or from Eq. (4) that 

H(k) « VA^[B. AB..... A^B] - VA*W 


( 9 ) 


where W is a controllability matrix (Ref. 18). Again if the system with order 
.! has no repeated eigenvalues, there exists a column vector B such that U has 
rank n. This means that H(k) is invertible if the system Is controllable and 
observable, letting k * 1. Eq. (8) will then determine the state matrix A In 
the following way 

VAV- 1 «H(1)H-1(0) (10) 

To rigorously prove this result, define E as the column vector E t *[ 1 ,0,...,0]. 
The measurement function y(k+l) can then be written by 

y(k+l) » ETH(k)E - ETH(k)H (0)H(0)E » ET[H(1)H (0)] k H(0)E (11) 

with the aid of Eqs. (8) and (10). Hence by Eq. (3), the triple [H(l)H“ l (0), 
H(0)E, ET] Is a realization in the sense that if the triple [A, B, C] in the 
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system equations (1) an'" (2) is replaced by the [H(1)H*1(0) , H(0)E, E T ], the 
measurement functions y(k) are reproduced as proved in Eq. (11). In other 
words, state variable equations (1) and (2) are transformed to the following 
equations 


x(k+l) « H(l)H-l(0)x(k) + H(0)Eu(k) 

(12) 

y(k) - E^k) 

(13) 

where 7(k) = V“!x(k). 

Let us summarize the case as follows. 

(14) 


A finite-dimensional, discrete time, linear time invariant dynamical system 
with a single input and a single output is realizable if the state matrix A has 
no repeated eigenvalues, and the system is controllable and observable. 

Hulti-input and Multi-output 

The time-domain description for this case is given by the pulse-response 
(or initial-state-reponse) function known as Markov parameter 

Y(k) * CA^B (or Y(k) - CA k [ Xl (0),x 2 (0) x m (0)]) (15) 

where xj(0) represents the ith set of initial condition and k is an integer. 
Note that B is a nxm matrix and C is a pxn matrix. The problem of system 
realization is that, given the functions Y(k), construct constant matrices [A, 
B, C] in terms of Y(k) such that the identities of Eq. (15) hold. The algorithm 
begins by forming the rxs block matrix (generalized Hankel matrix) 


H rs ( k— 1) 


rv(k) 

YUl+k) 


,Y(k+ti) 

.YUl+k+tx) 


Y(k+t s .i) 

, . . . ,Y( j^+k+tg_j ) 


|Y ( J r-1*^) (j r-l +lc+t l) • • • • (jr-l^^s-1 )J 


(16) 


where j-j (i *1 , . . . ,r-l ) and ti(i«l,...,s-l) are arbitrary integers. For the 
system with initial-state-response measurements, simply replace H rs (k-1) by 
H rs (k) % If Is easy to prove from Eq. (15) that Eq. (9) also holds for 
this multi-input and multi-output case. 


H rs (k) » ; V r 


CA 


and We - [B,A tr B A^B] 


CA Jr *l 


(17) 


where V r and W s are respectively the observability and controllability mai~ices 
in a general sense. Note that V r and W s are rectangular matrices with dimensions 
rp x n and n x ms respectively. Assume that there exist a matrix H* satisfying 
the relation 
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- I„ (14) 

where In is an identity matrix of order n. Define 0 p as a null matrix wlt'.i 

order p, E 1 = [Ip, 0 p ,..., 0 p ] and eJ » [I., fly...., Oy]. In view of Eqs. 

(16) and (18), the measurement function Y(k+1) can be obtained through either 
of two algorithms A1 and A2. Ihe algorithm A1 is 

Y(k*l) • EpH rs (k)E, • Ej* r » k H s H»V ( K s £, 

' E p T CY r *W s H*] k V r W s E,, 

• Ej[H rs (l)H*] k H,. s (0)^ 1 <’») 

and the algorithm A2 is 

Wl) - E^OOE, - 

■ EpV$C H,v r w sl k E. 

* EjH rs (0)[H»H rs (l)] k E. (20) 

fence, by Eq.(16), H„(0)E,, eJ] or [H*H«0). E». EjH rs (0)] 

is a realization. There fs no doubt that the matrix H* plays a major role in 
solutions (19) and (20). What Is H*? Observe that, from Eqs. (17) and (18), 

HrstOj^H^O) - - V r W s - H r$ (0) (21) 

H* Is a pseudo-inverse of H rs (0) in a general sense. When the rank of H rs (0) 

equals to the column number of H rs (0), then H**[[H rs (0)]^H rs (0)]“*[H rs (0)3*. 

If the rank equals to the row number, then H*<H rs (0)] T [H rs (0)[H rs (0)]*J"*. 
The matrix H rs (l)H* has been used In structural dynamics area to identify 
system modes and frequencies. 13 Both are special cases representing either 
single input or single output which can not realize a system that has repeated 
eigenvalues, or a noise-free system unless the system order Is a priori known. 
A general solution for H* is given below. 

For an nth order system, find the nonsingular matrices P and Q such 
that™. 11 

H fS (0) * PDQ T (22) 

where tho rpxn matrix P and the nxms matrix Q 1 are isometric matrices (all the 
columns are orthonormal), leaving the singular values of H rs (0) In the diagonal 
matrix D with positive elements [d^^f ••••dn]. The rank n of H rs (0) is 

determined by testing the singular values for zero (relative to desired 
accuracy) 1 ? which will be described In the next section. Define 

H rs (0)- PDQ T - [PD][Q T ] - P d Q T (23) 
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Each of the four matrices [pJ,Q*,Wc,vI] has rank and row number n. By 
Eq.(17) with k=0, 

v r w 5 = H rs (0) • P a q T (24) 

Multiplying on the left by P d and solving for Q* yields 

™s s < p d p d)’ lp d V r w s = Q T (25) 

T is nonsingular because if U * W S Q(Q T Q)”^* W $ Q, then TU* I by Eq. (25). Since 
TU = I = U T for nonsingular T ar.d *J then 

« s CQ(PjP d )-'pJ]V r ■ I„ (26) 


u eiKe» by Eq. (18) 

H # = [Q][(PjP d ) _1 Pj] = CQKtT 1 ? 1 ] * QP{ (27) 

The dimension of matrices Q and P^ with rank n are respectively msxn and 
nxrp. To this end, summarize the case as follows. 


A finite-dimensional , discrete time, linear time-invariant dynamical system 
with multi-input and multi-output is realizable in terms of the measurement 
function if the system is controllable and observable. 

Note that no restrictions on system eigenvalues are given for this case. 
In other words, this technique can realize a system with repeated eigenvalues. 
The system (1) with this realization will be transformed into the following 


equation 

x(k+l) * H rs (l)H*x(k) + H rs (0)E m u (28) 

y(k) * EjiT (k) (29) 

where x(k) = V s H*x(Y). Or (30) 

7(k+l) * H # H rs (l)7(k) + gi (31) 

y(k) * EpH rs (0)7(k) (32) 

where 7(k) = W s x(k) (33) 


The realizations (28)- (33) are not of minumum order, since the dimension of x 
is the number of either columns or rows of the matrix H rs (0) which is larger 
than the order n of the state matrix A for multi-input and multi-output cases. 

With the aid of Eqs.(17), (18) and (27), a minimum order of realization 
can be obtained from 
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Y(M) * EpH r s (k)E » * 

“ S p V i "'VV«r“s^ 

* E;H rs (0)«P5v r * k M s qpjH rs (0)t 1 , 

* EjH rs (0)Q[pjH rs (l)Q] k p5H rs (0)E m 

■ EpPdtP'nrst 1 >03 k Q T E„, (34) 

where Eq.(23) has been used to obtain the last equality. This is the basic 
formulation of realization for ERA. 

The triple [Pd H rs U)Q» Q T Em» E o p d^ * s a min ^ nw,n realization, since the 
order n of P*H rs (l)Q equals to the dimension of the state vector x. The same 
solution in a different form for the case where - tj > 1 (1«l,...,r-l) can 
be obtained by completely different approach as shown in Refs. [3 A 20]. The 
system (1) with this realization is written as 


x(k+l) * pjH rs (l)gx(k) + Q^u (35) 

y(k) - E T P d 7(k) (36) 

where 

x(k) = W s Qx(k) (37) 


A simple exercise such as replacing Y(k+1) by Y(k) in Eqs.(19), (20) and (34) 
shows that all the algorithms developed above are also true for the realization 
of a system with initial-state-response. 


Examination of Eqs. (19), (20) and (34) reveals that algorithms (Al) and 
(A2) are special cases of ERA. Al Is formulated by inserting the Identity 
matrix (18) into the right hand side of the state matrix A as shown In Eq. 
(19). On the other hand, A2 Is obtained by inserting the Identify matrix (18) 
into the left hand side of the state matrix A as shown in Eq. (20). However 
the algorithm ERA is formed by inserting the Identity matrix (18) Into both 
sides of the state matrix A as shown in Eq. (34). Because of the different 
insertion, Al and A2 do not minimize the order of the state transition matrix. 
Mathematically, if the singular val 3 decomposition technique is not included 
in the computational procedures, Al and A2 can not be numerically Implemented, 
unless a certain degree of artificial noise and/or system noise are present. 
Noises tend to make up the rank deficiency of the generalized Hankel matrix 
H rs (0) for algorithms Al and A2. Since the degree of noise presence is generally 
unknown, algorithrs Al and A2 are not recommended. 


III. MODAL PARAMETER IDENTIFICATION AND MODEL REDUCTION 


The presence of almost unavoidable noise and structural nonlinearity 
Introduces uncertainty about the rank of the generalized Hankel matrix and, 
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hence, about the dimension of resulting realization. By employing the singular 
value decomposition (SVD) technique, the rank structure of the Hankel matrix 
can be quantitatively displayed. Ihe set of singular values can be used to 
judge the distance of the matrix with determined order to a lower-order one. 
Therefore, the structure of the generalized Hankel matrix can be properly 
exploited to efficiently solve the realization problem. These include an 
excellent numerical performance, stability of the realization and flexibility 
in determining order-error tradeoff. 

Assume that, by Eq. (22) 

0® Diag. [dj, d2,...,d^,dp^i ,..., d|j] (38) 

with 

d H ^2 > ..* 1 d n >_ d n +l >,...>. dN (39) 

If the matrix H rs (0) has a rank n then all the singular values dj(i®n+l N) 

should be zero. When singular values di(i*n+l,...,N) are not exactly zero but 
very small, then one can easily recognize that the matrix H rs (0) is not far 
away from a n-rank matrix. However, there can be real difficulties in deter- 
mining a gap between the computed last nonzero singular value and what should 
be effectively considered zero, when measurement noise is present. Possible 
sources of the noise can be attributed to the measurement signal, computer 
round-off and Instrument imperfections. 

book at the singular value d n of the matrix H rs (0). Choose a number 6 
based on measurement errors Incurred in estimating the elements of H rs (0) and/or 
round-off errors incurred in a previous computation to get them. If 6 is 
chosen as "zero threshold" such that 6 < d n , then the matrix H rs (0) is 
considered to have rank n. Unless information about the certainty of the 
measurement data are given, the number 6 is defined as a function of the 
precision limit in the computer machine. For example, 6 *d n /di cannot exceed 
the precision limit, further details are found in Ref. [11]. 

After the test of sir >lar values, assume that the matrix [P^H rs (k)Q] has 
rank n. Find the eigenvalues c and eigenvectors i such that 

V’kPjHrsWQ^ a z (*0) 

The modal damping rates and damped natural frequencies are simply the real and 
imaginary parts of s, after transformation from the z- to the s- plane using 
the relationship 

s - C(lnz) + 2j* ]/(kAx) (41) 

where at is the data sampling interval and j Is an integer. Note that k is 
generally chosen as 1 for simplicity. Although z and are in complex domain, 
computation of Eq.(40) can be performed in the real domain (Ref. 21) sirce the 
state matrix realized for most flexible structures has Independent eigenvectors. 

The triple [ z, ♦” 1 Q T E JB , Ej[P d + ] Js obviously a minimum order of realiza- 
tion simaly by observing Eq. (3$). eIp^ Is called sensor modal displacements 
and +“'Q'E m initial modal amplitudes. To quantify the system modes and noise 
modes, two indicators are developed as follows. 
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Modal Amplitude Coherence y 


If the information about the uncertainties of the measurement Is minimum, 
the rank thus determined by the SVD becomes larger than the number of excited 
and observed system modes to represent the presence of noises In modal space. 
In modal parameter Identification, the Indicator referred to as modal amplitude 
coherence Is developed to quantitatively distinguish the system and noise 
modes. Based on the accuracy parameter, the degree of the modal excitation 
(controllability) Is estimated. 

Ihe modal amplitude coherence Is done by calculating the coherence between 
each modal amplitude history and an Idea one formed by extrapolating the Initial 
value of the history to latter points using the Identified eigenvalue, let 
the control input matrix (Initial condition) be expressed 

♦‘ 1 Q T E m - C^.bj, b n ]*, (42) 

where * means transpose and complex conjugate, and the 1 x m column vector bj 
corresponds to the system eigenvalue Sj( j«^,...,n) . Consider the sequence 

<ij = [bJ.expltjATSjJbJ exp(t s _ 1 AxSj)bj] (43) 

which represents the ideal modal amplitude In complex domain containing 
informations of the magnitude and phase angle with time step At . Now, define 
vectors qj such that 

♦ -1 Q T * Cqi»q2 q n 3* ( 44 ) 

The complex vector qj represents the modal amplitude time history from the real 
measurement data obtained by the decomposition of „ne Hankel matrix, let Vj 
be defined as the coherence parameter for the jth mode, satisfying the relation 

Y j * l^j<ljl/(lflj<ljllflj <, jl) 1/2 < 45 ) 

where | | represents the absolute_value. The parameter yj takes only the values 
between 0 and 1 . yj ♦ 1 as qj ♦ qj Indicates that the realized system 
eigenvalue Sj and the initial modal amplitude bj are very close to the true 
values for the jth mode of the system. On tne other hand. If yj Is far 
away from the value 1, the jth mode Is a noise mode. However, to make a clear 
cut between the system modes and noises requires further studies. Obviously, 
the parameter yj quantifies the degree to which the modes were excited by 
a specific Input, l.e. the degree of controllability. 

+ 

Modal Phase Collinearlty u 

lightly damped structure, normal mode behavior should be observed. An 
Indicator referred to as the modal phase collinearlty Is developed to measure 
the strength of linear functional relationship between the real part and the 
Imaginary part of the sensor modal displacement (mode shape) for each mode. 
Based on the accuracy Indicator, the degree of the modal observation Is 
estimated. Define 
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Ep^ ~ [Cj jCg , • . • »Cp3 (46) 

where Cj(j=l,2,...,n) is the sensor modal displacement corresponding to the 
jth realized mode, let the column vector 1 of order p be 

1 T » CM 1] (47) 

in which p is the number of sensors. Now compute the following quantities 
for the jth mode shape. 

Cj = cjl/p (48) 

c rr = [Rea1(Cj-Cjl)] T [Rea1(Cj-Cjl)] (49) 

c ri = [Real(Cj-Cjl)] T [Imag(Cj-cjl)] (50) 

c^ = [Imag(Cj-Cjl)] T [Imag(Cj-Cjl)] (51) 

e = (cn - c rr )/2c ri (52) 

© * arctan[e + sgn(e)(l + e 2 ) 1 / 2 ] (53) 

where Real( ) and Imag( ) respectively arc the real part and Imaginary part of 
the complex vector ( ), and sgn( ) is vne sign of the scalar ( ). The modal 
phase collinearity yj for the jth mode Is then defined as (Ref .22) 

Uj = f c rr + c ri [2(e 2 +l)sin 2 (e)-l]/e }/(c rr +c 1i ) ; j«l,2,...,n (54) 


This indicator checks the deviation from 0° - 180° behavior for components of 
jth identified sensor modal displacement. The parameter uj takes only the 
values between 0 and 1. pj ♦ 1 indicates that the accuracy of the modal 
displacement is high. On the other hand, if yj Is away from 1, the jth mode 
is either a noise mode or high damping Is present. 

Model Reduction 


The dynamical system is composed of an interconnection of all the ERA 
identified modes. The accuracy Indicators allow one to determine the degree of 
individual mode participation. Model reduction can then be made by truncating 
all the modes with low accuracy indicators. The accuracy of the complete modal 
decomposition process can be examined by comparing a reconstruction of Y(k) 
formed by Eq.(35) with the orginal free decay responses, using the reduced 
model . 


IV. SUMMARY OF ERA 


A flowchart of the procedures to be followed to use ERA In system model 
identification is presented in figure 1. The computational steps are summarized 
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as follows: 


1. Construct a block -Hankel matrix H rs (0) by arranging the measurement 

data into its rows with given r, s, t< (1*1,2 s-1) and (i * 1, 2 

r-1), (Eq. 16). 

2. Decompose H rs (0) using singular value decomposition (Eq. 23). 

3. Determine the order of the system by examining the singular values of 
the Hankel matrix H rs (0) (Eq. 38). 

4. Construct a minimum-order realization (A, B, C) using a shifted block- 
Hankel matrix (Eq. 34). 

5. Find eigensolutions of the realized state matrix (Eq. 40) and compute 
the modal damping rates and frequencies. (Eq. 41). 

6. Calculate the coherence parameter (Eq. 45) and the collinearlity 
parameter (Eq. 54) to quantify system modes and noise modes. 

7. Determine the reduced system model based on accuracy Indicators, 
reconstruct function Y(k) (Eq. 35) and compare with measurement data. 

Note that the determination of r, s, ti and In Step 1 above requires 
further development. This determination Is related to the choice of the 
measurement data to minimize the size of the Hankel matrix H rs (0) with the rank 
unchanged. 


V. EXAMPLES: SIMULATION AND EXPERIMENT RESULTS 


To illustrate the ERA method, two examples are given. First, a numerical 
problem will be posed and solved for an assuned structure with distinct and 
repeated frequencies. Second, experimental data for a simple, two-dimensional, 
grid structure as shown In Fig. 2 Is used and realized In terms of a linear 
system. Experimental results are compared with those predicted by a finite 
element model. 

Numerical Simulation 

Figure 2 shows a representation of a typical flexible structure. The 
dynamical equation for this typical structure with Initial-state-response in 


terms of system modes In modal space can be written as: 

dg/dt - A g (55) 

y ■ T g (56) 

where A is a canonical matrix with the diagonal blocks {A] A^}, 
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g is the generalized modal amplitude and C is the generalized sensor influence 
matrix. The quasi -diagonal matrix Aj (j*l, — ,k) has the matrix form 


A j 



(57) 


The complex values 6j + Iwj are the eigenvalue of the frame structure. 

Given a model described as in Eq. (55), results of some numerical simulation 
using the ERA scheme can be summarized in the sequel. Two cases will be given 
including systems with and without repeated eigenvalues. The numerical test is 
performed by taking as “data" y the output values of the solution of a model 
with the form (55) whose parameters A, C and initial condition g(t 0 ) are 
known. In the analysis of physical systems, experimental methods generate the 
measurement data y. It is then desired to realize a system by using the data y 
and compare with the known model. 


Case I: A model with distinct eigenvalues 

Assume that parameters such as bending rigidity, mass density and damping 
coefficient of the assumed structure are adjusted to give 


A j 8 


-O.Olxj j 

-j -O.Olxj 


j 


2, 3, 4, 5 


(58) 


To Illustrate applications of ERA in a single input and single output case, A 
sensor is chosen and located to give 

C = [l.O.i.O, 1,0,1»0,1,0] (59) 


Let the initial condition for free decay responses be 

g T (t 0 ) = [0,1,0. 1,0, 1.0, 1,0,1] (60) 

Then the functions y with a sample time Interval 0.05 second generated 
from the model (55) with known parameters (58), (59) and (60) are used as 
measurement data for the ERA procedure. 

Using = tj = i and r=s=90 in Eq. 16, the ERA realization of a dynamical 
system is 

C = [0.709,2.529,-0.347,-1.706,0.814,-1.183,-1.382,-0.276,1.129,1.257] (61) 

g T (t o )=[0. 103 ,0.367, -0.1 14, -0.563, 0.395,-0.574, -0.696, -0.139, 0.396, 0.440](62) 

and A is identical to that shown in £q. (58) with the accuracy close to the 
precision limit of the computer. In the process .af realization, the number 
6 *d n /d 1 as defined in Eq. (38) is set to be 10~'S The singular values of 
the generalized Hankel matrix H rs (0) are 
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D = [49.86,44.84,33.69,27 .64,23.69,21 .04,1 3.57,10.95,6.374.5.508] (63) 

All the values dj (j -11,..., 90) which has the number dj/d^ less than 10"^ 
are considered to be zero. The rank of the Hankel matrix J K rs (0) is obviously 
ten which is identical to the order a priori given in Eq.(58). The realized 
state matrix is a minimum order of 10 and the eigensolutions are obtained from 
this 10 x 10 matrix. All the parameters for modal amplitude coherence (Eq.45) 
and modal phase collinearity (Eq. 54) are 100%. Although Eqs. (61) and (62) 
are a different realization from the system (59) and (60), they are equivalent 
in the sense that a unitary transformation and normalization will make them 
equal . 

By forming the matrices V in Eq. (5) and W in Eq. (9) with the aids of 
Eqs.(58)-(60) , the reader can see that this realization is controllable and 
observable. 


Case II: A model with repeated eigenvalues and Independent eigenvectors 
Assume now that the system model is represented by 


and 


A l = a 2 


- 0.01 1 . 0 ' 
- 1.0 - 0.01 


<).01xj j 
-j -O.Olxj 


j* 3, 4, 5 


(64) 

(65) 


Using the same process as last case, the ERA realization simply miss the repeated 
eigenvalue K\. The result Is expected since, by control theory for a linear 
system, single input or single output does not make a system with repeated 
eigenvalues and independent eigenvectors controllable or observable. It can 
be verified that the matrices V in Eq. (5) and U in Eq. (9) formed by Eqs. 
(59), (60), (64) and (65) have rank 8. Multi-input and multi-output must be 
used to realize such a system. Let two sensors be chosen and located such 
that 

[ 1 , 0 , 1 , 0 , 1 , 0 , 1 , 0 , 1,01 


C = 


( 66 ) 


[ 0 , 0 , 1 , 0 , 1 , 0 , 1 , 0 , 1,0 


and two initial conditions for free decay responses 



" 0 , 1 , 0 , 1 , 0 , 1 , 0 , 1 , 0 , 1 " 
0 , 0 , 0 , 1 , 0 , 1 , 0 , 1 , 0 , 1 . 


(67) 


Note that the rows in Eqs. (66) and (67) are Independent. For each initial 
condition, a series of "measurement" function y with a sample time Interval 
0.05 second can be generated from the model (55) where each y in this case is a 
vector with two elements for two different sensors. The free decay function 
Y in Eq. (15) is then a 2 x 2 matrix. Using that ji * t^ * 1 and r * s * 45 
for Eq. (16), the ERA realization for a dynamical system is then 
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(68) 


C = 


0.135,-1.686,0.155,-0,172,0.111,-0.032,0.099,0.035,0.195,0.177 

-0.004,0.107,0.142,-0,136,0.111,-0.032,0.099,0.035,0.195,0.177 


g T (t 0 ) 


-0. 014, -0.457, 3, 840, -3. 692, 8. 338, -2. 406, 8. 9 F . ,3. 131, 2. 818, 2.554 
-0.051 , 0.092,3.605,-3.508,8.338,-2.406,8.956,3.181 ,2.818,2.554 


(69) 


where A is identical to that shown in Eqs. (64-65). The singular values Dare 
D = [70.16,44.32,37.97,25.25,1 1.18,9.050.7.950,3.873,0.1 27,0.026] (70) 


The seme error window 6=d n /d-j as last is used. All the parameters for modal 
amplitude coherence and modal phase colinearity are 100%, Again, Eqs. (68-69) 
and Eqs. (66-67) are equivalent in the sense that a unitary transformation and 
normalization will make them equal. The reader can easily verify that this 
realization is controllable and observable. 


Sample Experimental Results 

A sample set of modal Identification results that have been obtained from 
laboratory .est data using ERA are included in this section. The test article, 
shown in Fig. 2, is a 7 ft by 10 ft flexible grid structure that will be used at 
NASA Langley for vibration control experimentation. It is constructed of over- 
lapping aluminum bars of 1/4 In. by 2 in. cross section, riveted together at 
one-foot intervals. Four rivets are used at each joint to provide a tight 
connection. The structure is sc -pended from a stiff overhead beam using two 
short cables attached to the top horizontal member. The results to be shown 
are from a preliminary dynamics test of the grid. It was conducted by exciting 
the structure with an air jet and measuring the free vibration response using 
nine non-contacting proximity sensors. The response sensors were attached to 
a stiff frame located adjacent to the grid for the measurement of out-of-plane 
motions. Eight different excitation frequencies corresponding to resonant 
responses were used. The sampling rate was 32 samples per second. 

The ERA analysis was performed using a single matrix of data from all 
nine response measurements and eight Initial conditions. Each response func- 
tion Y as shown in Eq. (16) was thus a 9 x 8 matrix. The Hankel matrix H rs 
of 72 rows by 400 columns was formed to perform the analysis. Table 1 provides 
h comparison rf the ERA results with analytical prediction from a NASTRAN 
finite-element model. The entries in the center of table are correlation 
coefficients in percent between each ERA-1 dentlfied mode shape and each NASTRAN 
mode shape. High correlation values Indicate good agreement between the two 
shapes. The results show reasonable agreement In both frequencies and mode 
shapes, except for the damping result of the first mode. The main reason for 
the first mode discrepancy is inadequate data length. Only 50 data points were 
used which corresponds to less than one cycle of data for the first mode. The 
results can be improved by using more data points. Note that few high correla- 
tions occur for some modes with significantly different frequencies. This Is 
because only 9 sensors were used In comparison. More detailed experimental 
results for a complex structure are shewn in Ref. [19], 
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CONCLUDING REMARKS 


An Eigensystem Realization Algorithm (ERA) is developed for parameter 
identification and model reduction for dynamical systems. Twi developments 
are given in this paper. First, a new approach is developed to derive the 
basic ERA formulation of minimum realization for dynamical systems. As by- 
products of this approach, two alternative less powerful algorithms, identi- 
fied as A1 and A2, are derived. A special case of A1 Is shown to be equiva- 
lent to an approach currently in use in structural dynamics. Second, accu- 
racy indicators are developed to quantify the partlpatlon of system ■"C^es 
and noise modes in the realized system model. In other words, degree of 
controllability and observability for each participated mode is determined. 
A model reduction can then be made for controller design. 

important features of the ERA algorithm are summarized as follows. 

(1) From the computational standpoint, the algorithm is attractive, since 
only simple numerical operations are needed; (2) the computational procedure 
is numerically stable; (3) the structural dynamics requirements for modal 
parameter identification and the control design requirements for a reduced 
state space model are satisfied; (4) data from more than one test can be used 
simultaneously to efficiently identify the closely spaced eigenvalues; (5) no 
restrictions on number of measurements are imposed. 
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Table 1 : Comparison of the ERA results with the NASTRAN Model 
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Free-response functions 
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Figure 1, Flow Chart of FRA 
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